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Abstract. The dissipative quantum dynamics of an anharmonic oscillator is 
investigated theoretically in the context of carbon-based nano-mechanical systems. 
In the short-time limit, it is known that macroscopic superposition states appear for 
such oscillators. In the long-time limit, single and multi-phonon dissipation lead to 
decoherence of the non-classical states. However, at zero temperature, as a result 
of two-phonon losses the quantum oscillator eventually evolves into a non-classical 
stationary state. The relaxation of this state due to thermal excitations and one- 
phonon losses is numerically and analytically studied. The possibility of verifying the 
occurrence of the non-classical state is investigated and signatures of the quantum 
features arising in a ring-down setup are presented. The feasibility of the verification 
scheme is discussed in the context of quantum nano-mechanical systems. 
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1. Introduction 

Nanoelectro-mechanical (NEM) resonators such as cantilevers, membranes and beams 
are important systems for the study of quantum phenomena in macroscopic, mechanical 
man-made objects |1|. For about a decade, the goal has been to cool the NEM resonators 
to the ground state and to verify the achievement of this state. This has recently 
been accomplished (2]-[4], and the next step is the generation, detection and coherent 
manipulation of more intricate quantum states, e.g., superpositions of macroscopic states 
(Schrodinger cat states) |5]. 

While a harmonic quantum oscillator displays a behavior analogous to its classical 
counterpart, the same does not hold for nonlinear oscillators. The presence of a 
conservative (Kerr) nonlinearity facilitates the creation of non-classical states. Although 
top-down fabricated NEM-resonators typically possess a finite Kerr constant /x, it is 
usually too small to alone suffice for non-classical state generation. In contrast, carbon 
nanotubes (CNTs) and graphene membranes can exhibit strong enough nonlinearities for 
non-classical states to emerge during free evolution (5[|6j. Using CNTs or graphene also 
alleviates the need for active cooling due to their low masses and high frequencies |1,5|. 
For this reason, carbon based NEM-resonators are highly promising for the study of 
nonlinear mechanical systems in the quantum regime. 

If a system possesses conservative nonlinearities, it is natural to also take into 
account higher orders in the system-environment coupling. This results in nonlinear 
damping (NLD) (T]- 10 . In carbon based NEM-resonators it has recently been shown that 
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the dissipative mechanisms can indeed be nonlinear in the oscillator coordinate 
Hence, to fully understand the quantum dynamics of carbon based NEM resonators it is 
essential to comprehend the interplay between conservative nonlinearities and NLD. In 



general, NLD will create interesting quantum features rather than remove them 13 -17 



For example, as is well known in quantum optics [13-15, 18 , for a harmonic oscillator, 
the stationary state obtained from NLD due to two-photon losses at zero temperature 
(T = 0) is a mixed state. Moreover, because of the parity conserving nature of the 
NLD, this stationary density matrix p has coherences in terms of nonzero off-diagonal 
elements poi in the number basis. Recently, it was proposed to use this effect in order 
to generate and preserve Schrodinger cat states in a double well potential realized by a 
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In this paper we study a situation more relevant to the coherent quantum dynamics 
of carbon based NEM-resonators, namely, a nonlinear oscillator with NLD. In this field 
there is little to be found in the literature |9|. Our main goals are to characterize 
the resulting stationary states along with the corresponding features of the relaxation 
dynamics. For the nonlinear oscillator with NLD we find that the parity conserving 
nature of NLD helps to preserve phase coherence also in the presence of a finite 
conservative nonlinearity. However, we also find that for sufficiently large ratio between 
the Kerr constant /x and the NLD strength 72, the coherence of the stationary state is 
suppressed. Further, at finite temperatures, with NLD, the final state is not a standard 
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thermal state with a reduced density matrix p oc e~^/^^^ unless linear damping (LD) 
is present. The corresponding Wigner distribution has a dip in the center, whose depth 
depends on the initial conditions. In the presence of both NLD and LD, the relaxation 
is governed by a crossover between two different relaxation rates. In the long-time limit, 
the decay rate V is independent of the initial oscillator amplitude for NLD. Instead 
it is a function only of temperature T. To facilitate experimental determination of the 
governing damping mechanism, we also present the corresponding signatures in the time 
evolution of the position coordinate variance. 

This paper is organized as follows: In section [2} we present the underlying 
theoretical frameworks used for analytical (rotating wave approximation [RWA]) and 



numerical (Wangsness-Bloch-Redfield 19 -■21]) calculations. The latter are used to verify 
the validity of the RWA-calculations. Then, in section [3[ we present results pertaining 
to the stationary states for various scenarios involving NLD and LD. Results on the 
relaxation towards the stationary states are found in section [4| followed by a section 
with summary and conclusions. 

2. Model and Methods 

2.1. Model 

To model the dynamics of the anharmonic oscillator subject to linear and nonlinear 
damping, we consider a single mechanical mode with frequency Vt and Kerr constant /x. 
The mode is coupled linearly and nonlinearly to a reservoir of harmonic oscillators with 
mode frequencies ujk- The mode mass and the reduced Planck constant are set to be 
m = 1 and h = 1^ respectively. Correspondingly, the total Hamiltonian is 

H =Hs + Hb + Hsb, (la) 

Hs =^/ + i^2V + fg^ (lb) 

= ^oJkblbk, (Ic) 
k 

HsB = g J](^77i,)(6l + bk) + J2(^rj2k){bi + b^), (Id) 

k k 

where the mode position and momentum operators are q = qo{a^ + a)/V^, p = 
i{a'^ — a)/v^go5 similarly the bath operators are given by Xk = ^o/c(&I + bk)/V^^ 
Pk = i{bl — bk)/V2xok- Here qo = and Xok = ^/^/'i^k^k are the zero point 

amplitudes of the oscillator and the bath, and satisfy the canonical commutation 
relations. The linear and the nonlinear coupling constants of the fc'th reservoir mode to 
the resonator mode are denoted by rjik and r]2k^ respectively. These couplings lead to 
amplitude-independent (linear) and amplitude-dependent (non-linear) [7-10 damping 



of the oscillator motion. The former coupling involves the exchange of single vibration 
quanta (vibrons) with the reservoir, while in the latter case two vibrons are transferred 
simultaneously. 
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We consider the following scenario: The oscillator mode is assumed to be initially 
in the ground state |0). Then, by applying a voltage pulse for a finite time, the 
oscillator is displaced. The corresponding state is, by definition ^22j, a coherent state, 
^{t = 0) = D[a) |0) = \a). Here, D{a) = exp{aa^ — a^a) is the displacement operator 
and the amplitude a depends on the strength and duration of the voltage pulse. During 
dissipation-free evolution under only, a coherent state will evolve into a Yurke-Stoler 
cat state {\a)+i\—a))/ \f2. This state will periodically re- appear with a frequency set by 
the Kerr constant /x. Typically the state's emergence is slow compared to the system's 
eigenfrequency Vt |!23|. When taking into account the interaction with the reservoir, the 
created cat states will be infiuenced by the dissipation and the oscillator will eventually 
relax into a stationary state. The state of the mechanical mode during the relaxation 
is described by the reduced density matrix (RDM) p = tiB Pt, which is obtained from 
the state of the total system px by tracing over bath degrees of freedom. 

The time-evolution of the reduced density matrix is described by a (generalized) 
quantum master equation (QME) [24] . Here we follow the standard approach using the 
Born-Markov approximation to obtain the QME from the von Neumann equation for 



the total density matrix 24 . In this approximation the QME in the interaction picture 
with respect to Hs is given by 

d r 1 

= [Mt)Amit - r)p{t) - A^{t - T)p{t)Ai{t)\ a^(r) 

+ [ p{t)A^{t - T)Ai{t) - Ai{t)p{t)A^{t - r)] C^i{-t). (2) 

The operators A^ are given by the expansion of the interaction Hamihonian in the 
interaction picture 

HsB = Yl ^-(^) ® ^-(^) ' (3) 

m=l,2 

where 

Am{t) = e'""'' {a^+ar e-^^^^ B^t) = J] q^r]mk {ble'^'' + bj,e-'-'') .(4) 

k 

Assuming the reservoir to be initially in thermal equilibrium, the reservoir correlation 
functions are given by 

Cmi{r) = tTB{Bm{t)Bi{t-r)pB} 



/I 
^kU^) K(u;)e''^" + {n^{u) + 1)6-^^] , 



(5) 



where nB{u)) = [e^/^^^ — 1] ^ is the Bose distribution function, pe is the thermal state 
and K.rni{^) = i^im{^) = "^TT Y,kiVmkVikQ^^^)^{^ - ^k) IS the spcctral density of the 
environment coupling. 

In general, the frequency dependence of tv^i is determined by the specific geometry 
of the system. As long as the spectral density is smooth and slowly varying around 
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UJ = and uj = 20. the detailed dependence on uj is not important [2^ (see also 2.2). 
To be specific, we consider an Ohmic spectral density |24|, i.e., we take 

= ^ml ^r—- . (6) 

iWml 



The values of Fn and are used to specify the strength of LD and NLD, respectively. 
The off-diagonal value can, in general, be chosen independently of Fn and r22. If the 
oscillator couples to two independent baths, the off-diagonal spectral density vanishes. 
In the RWA, which is discussed in |2.2| the dissipation terms associated with ri2 can 
be neglected. Here we set ri2 = a/TiiT^, which vanishes if either linear or nonlinear 
damping is absent. 

2.2. Relaxation dynamics in the rotating frame 

Typically, for NEMS the time scales associated with Q and fi are well separated with 
/X (see also section [5]), which justifies an analytic treatment in the RWA. Within 
this approximation the QME ^ in the Schrodinger representation becomes |26| 

dtp{t) = {Co + Ci + C2)p{t), (7) 

where the Lindblad superoperators are given by 

CqP = —i{Q + iii)[a'^a^p]—iiii[{a'^ay^p]^ 

^iP = ^-apa'^ — ^{a'^a, p} + ^^a'^ pa — ^{aa'^ , p}, (8) 
^2P = — ^{a'^a^aa^ p} + ct^ pctct — ^{aaa'^a'^ ^ p}. 

The superoperators Cq^ C\ and C2 describe the coherent evolution due to i/s, linear 
and the nonlinear interaction with the reservoir, respectively. The linear and nonlinear 
dissipation rates, 7± and 72^, are Fourier transforms of the correlation function ([5]), and 
are given by 

7+ = riinB(f^), 72+ = r22nB(2f]), 

7- = rn[nB(f^) + 1], 72- = ^22\n^{m + !]• 
We would like to note, that the RWA equations given above are only valid for a weak 



nonlinearity 26 . This is reflected by the dependence of the rates on the transition 
frequency of the harmonic part of the Hamiltonian i^s and not on the frequencies of 
the full nonlinear i/s- The validity of this approximation in the present case will be 
shown in the next sections by comparing to numerical calculations which include the 



exact transition frequencies (see also section 2.3). 

The solution to ^ is obtained by rewriting the matrix equation as a set of decoupled 
equations. To this end, p is projected onto the number basis Pn,n+\ (n|p|n + A), where 
A = 0, 1, 2, . . .. It is convenient to define pn,n+\{t) = ^J^^'^^ni^ t)e^^(^+^)* 13 , which 
also removes oscillations due to the term linear in a^a in ([t]). For fixed A the function 
'0^(A,t) gives the entries of the A'th super diagonal of p. Hence, ([T]) becomes 

dtipni^^t) = 72+(n + A)(n + A - l)V^n-2 + 72-(^ + 2)(n + l)V^^+2 
- {-i/xA(A + 2n) + 72-n(n - 1 + A) 
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+ ^A(A - 1) + 72+(n + l)(n + 2 + A) + ^A(A + I)}./;, 
+ 7_(n + + 7+(^ + 

- 7-(^ + ^A)^/;, - 7+(n + 1 + ^A)^/;,. (10) 

As indicated above, the equations for different A do not couple. In particular, for A = 0, 
we obtain a single equation for the populations Pn,n{t) = '0n(O,t), which is used in the 
following to analyze the relaxation to the stationary state. 



In general, when all coefficients in (10) are non-zero the time evolution of '0^(A,t) 



has to be found numerically. If either 7_ or 72- is vanishing, i.e., there are only one or 
two vibron losses present, the time-dependence of ipn{^i t) can be calculated analytically 



(see 27-29 and [13]). In these cases it is also straight-forward to find the stationary 
solutions for all A. At finite temperatures, when in addition to the losses, one or two 
vibron excitations (7^ > or 72^ > 0) contribute, the stationary probability distribution 



Pn^n can be found using the detailed balance condition [30[[3l] (see Appendix B). If one 
and two- vibron processes compete (7^ > and 72^ > 0), a solution for p^^^ can be given 
in terms of a continued fraction |32| or (for 72+ = 0) in terms of confluent hypergeometric 
functions [33 1. 

In the next sections we are in particular interested in the behavior of the position 
variance 

<Ag^> = <g^> - {qf , (11) 

where (•) = trs[*p]. Expressing q in terms of creation and annihilation operators in the 
rotating frame one finds 

<Ag2> ^ I [<ata + aa^) - 2| <at> = 1 + {a^a) - \ <at> p , (12) 

where all oscillating contributions have been discarded. The expectation values of a'^a 
and a'^ can be expressed in terms of '0^(O,t) and '0^(l,t), respectively, 

oo oo 

(a^a) = "^npn^n = J]]n?/;^(0,t) , (13) 

oo oo 

(at) = Vn + lpn,n+l = V^n(l, t) . (14) 

n=0 n=0 

Knowing '0ri(O, t) and ^0^(1, t) therefore allows for an evaluation of the position variance. 
In the following we focus on obtaining such solutions in the long-time limit. 

2.3. Numerical Computation 

In order to investigate the oscillator dynamics and to verify the results obtained by 
the RWA approximation, the QME ^ is solved numerically by Wangsness-Bloch- 



Redfield 19-21 calculations. To this end the operators are represented in the 
eigenbasis of the oscillator Hamiltonian Hs and the Hilbert space is truncated after 30 
states (40 states for a = 3), which yields converged results. The initial coherent state. 
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p(0) = |o^)(Qf|, is generated by applying the displacement operator D{a) to the ground 
state. 

3. Stationary state 

3.1. Zero temperature 

We begin with the RWA-analysis of the stationary state p(oc) at zero temperature. To 
evaluate the applicability of the RWA, we compare the results to numerical calculations 



as described in section 2.3, If only single vibron losses are present in the system (7+ = 0, 
l2± = 0), one finds from the equation of motion ([To]) for p that any initial state p(0) 
will decay to the ground state, i.e., p(oc) = |0)(0|. The associated position variance is 
(Ag2) /qi = 1/2. 

On the other hand, if only nonlinear losses are present (72+ = 0, 7± = 0), a general 
initial state will relax to the stationary state fl3||l4] 

pM=Po,o|0)(0|+Po,i|0)(1|+Pi,o|1)(0|+Pi,i|1)(1|, (15) 

where the matrix elements in ( [T5| ) depend on the initial density matrix p(0). Since for 
two- vibron losses parity is conserved, the only nonzero matrix elements on the main 
diagonal of p(oo) are 

PO,0(OC) = V^o(0, 00)=^ Pn,n(0) = Peven, (16) 

n even 

Pl,l(00) = ^1(0, 00)= ^ Pn,n(0) = Podd • 

n odd 

The stationary off-diagonal elements po,i and pi^o are in general non-zero and can also 



be found from a conservation law. From (10) one gets 



PO,l(oo) = pt^o(oo) = X]C2„^2n(l,0)e'("+2'^)* = 



n=0 



where the coefRcients C2n are given in terms of Gamma functions (see Appendix A) 



r(| + n)r(i-i:^) 

This expression is a generalization of the result of Simaan and Loudon |13| to the case 
of an anharmonic oscillator (/x 7^ 0). Note, that unless /x = 0, the ofl[-diagonal element 
Po,i is time-dependent in the stationary limit. Specifically, for a system which is initially 
in a coherent state \ip{t = 0)) = i.e., with matrix elements 



Pn,m(0) = exp(-|a|2)a^a*^/y^I^ , (18) 
the constants of motion are 

Peven = exp(- |a|^) COSh( |a|^) , 

^odd = exp(-|a|^) sinh(|a|^) , 
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Figure 1. Wigner distributions (20) of the two lowest number states and stationary 
states resulting from nonlinear ly damped, initial coherent states with a = 1 at T = 0. 
(a) Ground state |0)(0|. The black line indicates the size of the state given by 
^/{Aq^ /qo = 1/4. (b) First excited number state |1)(1|. The white line indicates 
the size of the ground state, (c) Weak NLD (72- //i = 1/10). The negative domains in 
the Wigner distribution reflect the presence of coherences (non-zero off-diagonals) in 
the stationary state. For comparison the black line shows the size of the ground state, 
(d) Strong NLD (72- //i = 10). The black line is included for the same purpose as in 
(c). 



and 

'0even(i) = Q^* exp(-|a|^ - i/xt)oFi(l - i— , |q^|^/4) . (19) 

72- 

Here, 0^1 is the confluent hypergeometric limit function [34]. For /i = the latter 
becomes a Bessel function of the flrst kind and the expression given in 13 is recovered. 
To illustrate the differences between the numerically calculated states obtained by 
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Figure 2. Zero temperature position variance of the stationary state resulting from 
NLD as a function of the displacement amplitude a of the initial coherent state. The 
variance is averaged over one period 2i\ j \i. Dashed lines show the behavior according 



to (21) for different values of the ratio 72- //i. Symbols denote the corresponding 
numerical results. The lower dotted line (red) indicates the variance of the ground 
state, resulting from LD. The upper dotted line (gray) indicates the upper variance 
limit, resulting from weak NLD for a ^ 00. 



LD and NLD, their Wigner distributions [35] 

W{q,v) = (2^)"'/ e-'^«(c? + e/2|p(oo)|(7-e/2), (20) 

are displayed in figure [TJ For comparison also the Wigner distributions of the ground and 
the first excited number state are included in figures [l] (a) and (b). While the ground 
state Wigner distribution is a Gaussian minimum uncertainty state, the first excited 
number state is a quantum state with negative domains in the Wigner distribution. 
The p(oc) resulting from NLD is a non-classical state with both negative and positive 
values in the Wigner distribution. This state is a mixed state which possesses coherences 
in terms of nonzero off-diagonal elements po,i in the number basis (cf. 
(c) and (d) display this stationary state for two different values of the ratio between the 
nonlinear damping rate and the Kerr constant. By comparing the sizes of the stationary 
states obtained by LD and NLD, it is seen that the latter occupy a larger phase-space 
area. This suggests, that the position variance (Ag^) is suitable for quantifying the 
differences between the states resulting from these damping mechanisms. 

Using the definition ( [T2| ) and the results above, one can evaluate the variance to be 

{Aq\a))T=o/q! = ^+e-l"l%inh(|an-|e-l"l'a*oFi(l-W72-, |a|V4)|'.(21) 

In figure |2| the variance according to ( [2T| ) is shown as function of the displacement 
amplitude a of initial coherent states for different values of 72- (dashed lines). The 



15)). Figures [T] 
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Figure 3. (a) Wigner distribution (20) of the stationary state resulting from NLD 
at finite temperature with parameters 72- = /i, a = 1, k^T/Vt = 1. A sector of the 
Wigner distribution is cut out for a better visuahsation of the dip in the middle, which 
indicates the difference to a regular thermal state, (b) Probability distribution of the 
even (dark blue bars) and odd (light blue bars) number state sectors of the stationary 
state shown in (a). The narrow bars represent the Boltzmann distribution of a thermal 



state described by (22). The insets show the Wigner distributions of the even and odd 
sectors, respectively. 



numerical results (symbols) are in accordance with (21). For comparison, also the 
variance of the ground state is included (red line), indicating its independence on the 
initial displacement amplitude. For small initial displacements, (a ^ 1), the behavior 
of the variance is governed by the coherences po,i = Pio^ since pi^i 0. For large 
amplitudes, (a 1), and weak NLD, (72- <C /x), the variance saturates at 1 (gray 
dotted line). This is due to the ground and the excited states being equally populated, 
Po,o = = 1/25 and the vanishing off-diagonal elements, since the hypergeometric 
function converges to zero. In the limit of strong NLD, (72- ^ /x), and large amplitudes, 
(a 1), the hypergeometric function converges to a non-zero value and the variance 
becomes (Aq^) ^ 1 — l/27r, which is larger than the variance of the ground state. 

3.2. Finite temperature 

Also at finite temperature p(oc) obtained by NLD will differ from the state obtained 
by LD. For LD only, p(oc) is fully characterized by the values of the diagonal matrix 
elements of the density matrix, 

p^,^ = Z{Q)-\-^''/^'-^\ (22) 

where Z{Q) = '^^exp[—nQ/{kBT)]. Here, /cb is the Boltzmann constant and T is the 



temperature of the bath. This result is found from (10) by taking A = and setting 



the time- derivative to zero. One obtains the detailed balance condition, which leads 



Multi-phonon relaxation and generation of quantum states in a nonlinear mechanical oscillatorll 



to (22). As we have pointed out in section 2.2, the RWA equations ([T]) are vahd for 
weak nonhnearities since the bath is only probed at the frequency Consequently, at 
finite temperatures the stationary probability distribution calculated within the RWA 
will not depend on the Kerr constant. However, for sufficiently small temperatures, i.e., 
/i (a'^'a) <C + /X, the approximation is valid as we will also show in the following. 

Since for nonlinear damping the two-vibron exchange conserves parity, one finds 
that detailed balance holds for each parity sector individually in this case. It can be 
shown (Appendix B) that 

P2n, 2n = Z (2^)-' e-^^""/ ^''^^^ P,,,^ , (23) 

P2n+l,2n+l = Z {2n)-' 6-'^''/^''^^^ P,dd ■ 

Note that the stationary probability distribution at finite temperature also depends on 
the initial state. This implies that expectation values will in general also depend on 
Peven and Podd- For example, the average number of vibrons in the oscillator mode is 
given by 

oo 

(at a) = npn^n = 2nB(2J^) + Podd , (24) 

n=0 



which is derived in [Appendix B 



Typical features of a stationary state obtained from nonlinear damping at finite 
temperature are displayed in figure [3j Figure |3] (a) shows the simulated Wigner 
distribution (see (20)), with parameters 72- = /x, a = 1, k^T/Vt = 1. A sector is 
cut out for better visualisation of its features. As discussed in section [3^ the thermal 
probability distribution of a nonlinearly damped state is split into an odd and an even 
parity sector. Each of the sectors has a Boltzmann-like distribution (see ([23])). The 
shape of the Wigner distribution refiects this parity dependence by having a dip centered 
at the origin. It can be shown that Vl^(0, 0) = (Peven — ^odd)/27r. Therefore, the dip 
in the Wigner distribution only depends on the initial state and not on temperature. 
Hence, in the case of an initial coherent state the dip will only depend on the initial 
oscillator displacement. For large displacement amplitudes the difference between the 
sectors becomes smaller, and Vl^(0, 0) 0. If the initial state had an odd parity, the 
dip of the Wigner distribution would attain a negative value. 

Figure [3] (b) shows the Boltzmann-like distributions of the odd and the even sectors 
of the stationary state discussed above. The narrow bars represent the distribution of 
the thermal state obtained by linear damping (see (22)). The insets show the Wigner 
distributions of the even and odd sectors. 

As the Wigner distribution of the finite temperature stationary state displays 
quantum signatures, so does its variance. Figure [4] shows the stationary state variance 
as function of the displacement amplitude a of the initial coherent state for several 
temperatures {k^T/Vt G [0,1/4,1/2,3/4,1] and 72- = /x). Dashed lines correspond to 
the analytical expression given by (33). Numerical calculations (symbols) are shown 
for the same parameters. The variance is seen to always be larger than the zero 
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Figure 4. Finite temperature position variance of the stationary state as a function 
of the displacement amphtude a of the initial coherent state. The variance is averaged 
over one period 2i\ j \i. Dashed lines show the behavior according to (33) for 72- = \i 



and different temperatures. Symbols denote the corresponding numerical results. 



temperature variance, displayed in figure [2} This is due to an additional contribution 
from the Boltzmann distribution and the vanishing off-diagonal elements of the density 
matrix (see section 14. 3|) . The variance is still dependent on the initial displacement and 



saturates at {/S.q^) /q^ = 1 + > 1 for large a, which is in contrast to k^T/Vt = for 
which {Aq^)/ql < 1. 

4. Relaxation Dynamics 

4.1. Zero temperature^ NLD 

We now investigate the relaxation dynamics of initial coherent states which evolve 
into the stationary states discussed in the previous section. Figure [5] (a) shows the 
numerically calculated variance (Ag^) as function of rescaled time r = f]t/(27r) for 
several values of the ratio 72-//^. The variance is sampled at multiples of the period 
f]/(27r), and therefore fast oscillations are not seen. For weak NLD (72_//x < 10), the 
short time dynamics {t < n/fi = 500) is governed by the Kerr constant /i. This is 
reflected in the initial increase of variance. The maximum variance, at which the Yurke- 
Stoler cat state is first expected to appear ^^||, is attained at r = 7r/(2/x) ^ 250. 
The variance of this state is given by (Ag^)cat = Qoi^^^ + l)/2, which is always larger 
than the variance of a coherent state for a nonzero a. As seen in the figure [s] (a), the 
shape of the variance peaks changes in time as the cat states decohere by NLD. Figure [5] 
(b) shows the Wigner distribution sampled at r = 7r/(2/x), displaying a cat state under 
infiuence of NLD. Some of its original features are still visible, like the remnants of the 
negative interference domains. 



After r = 7r//i, the system gradually relaxes to the state (15). In figure (a) this is 
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Figure 5. (a) Time evolution of the position variance of an initial coherent state 
with a = 1 at T = for different NLD strengths. The dashed line (red) indicates 
the variance of the ground state, (b) Wigner distribution W{q^p) of the decohered 
Yurke-Stoler cat state attained at r = 7r/(2/i). (c) Wigner distribution W{q^p) of the 
state in (b) at a later time, r = 37r//i, reaching the stationary state. Black contours 
encircle domains where < 0. 



reflected by steady oscillations around a mean value of (Ag^) which is, as discussed in 
section 3.1, larger than 1/2. The oscillations in (Ag^) are caused by the complex valued, 



oscillating coherences of po,i- The relaxation dynamics described above becomes faster 
with increasing ratio 72- /m- Figure [5] (c) shows the Wigner distribution of the state 
sampled at r = 37r/(/x), where the resemblance to the stationary state in figure [l] (c) is 
clearly seen. 



4^2. Zero temperature, LD and NLD 

As shown in the previous section, NLD at T = leads to the formation of the non- 



classical state (15). In the additional presence of LD the generation of this state will 
be affected and eventually the ground state will be reached. In the following we discuss 
the resulting decay of the variance, which can be used to extract information about the 
nature of the damping mechanism and the state initially created by NLD. We introduce 
the ratio 7_/72_ as a measure of the relative dissipation strength, and set 72- = /x. 



In presence of LD at T = 0, (10) leads to (for A = 0, 1 



d 

Q^PnAt) =l-{n+ l)Pn+l,n+l + 72- + 2)(n + l)Pn+2,n+2 

- [72-^(^ - 1) + 7-^] Pn,n , (25) 

d 

Q^^^ni^^ t) = 7-(n + l)V^n+i + 72- (^ + 2)(n + l)V^n+2 
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Figure 6. (a) Time-dependence of the position variance of a linearly and nonlinear ly 
damped mechanical mode for a = 1 and 72- = /i. Full lines denote different 
values of 7_/72_. Dashed lines show the function (27). (b) Time-dependence of 
(1 — po,o) with parameters corresponding to those of the main figure, (c) Time 
dependence of {Aq'^{r))/qQ — 1/2 for 72- = /i and 7_/72_ = 1/10 for different initial 
displacement amplitudes. Dashed lines indicate the behavior according to (27). For 
visual convenience, the graphs corresponding to a = 3 are multiplied by a factor of 2. 



-i/x(l + 2n) + + j-{n + ^) 



. (26) 



It is seen that only po,o is non-zero in the stationary limit, i.e., pi^i and ipn decay 
for sufficiently large times. In the limit of strong NLD, 7_ <C 72-, one finds that 
the populations p^^n for n > 1 rapidly decay due to two-vibron losses, but Pi is 
only decreased via one-vibron losses. Hence, pi^i ^ — 7_pi^i and ^ e~^-^Podd- 

Here it is assumed that the initial decay on the time-scale I/72- is not infiuenced by 
the LD and po,o(0) ^ ^even and Pi,i(0) ^ Podd- Similarly, one finds that '0o(l,t) ^ 
exp[(i/x — 7_/2)t]'0even(O). Therefore, the position variance will decay to the minimum 
uncertainty value 1/2 according to 

{Aq') /ql = i + [Podd - |^eve„(0)n e-^-* . (27) 

Figure [g] (a) shows the time evolution of the numerically calculated (Ag^) for 
different values of the ratio 7_/72_ and ji = 72-. For 7_/72_ ^ 1/10, the evolution of 
the variance for short times t < 7r//i resembles the result shown in figure [5j The short- 
time dynamics is again infiuenced by the Kerr constant, while the long-time behavior, 
(r ^ vr//x) is showing the relaxation of the state created by NLD. In this stage the 
variance is exponentially decaying with the decay being faster for larger values of 7-/72- . 
In figure |6] (b) the time evolution of (1 — po,o) is shown, confirming the decay to the 
ground state. The parameters are the same as in|6](a). For equally strong LD and NLD, 
(7_ = 72-), an exponential decay is seen for the entire time interval. For LD weaker 
than NLD, (7_ < 72-) and r < 7r/2/x, the decay is initially not exponential, whereas for 
r ^ n/iLi the ground state population displays exponential increase. 

Figure [g] (a) suggests the concept of sequential relaxation for weak LD. The first 
stage of the sequence (0 < r < 7r//x), is mainly dominated by NLD and the Kerr 
constant, which leads to the formation of the state in (15). The second stage (r > 7r//x), 
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is mainly influenced by the relaxation to the ground state because of LD. In other words, 



the state in (15) is initially created by NLD, and then decays due to LD. Accordingly, the 



average variance in the second relaxation stage is expected to be given by (27), which 



is displayed in flgure |6] (a) (dashed lines). It can be seen that the idea of sequential 
relaxation is quantitatively correct for 7_ <C 72-. Furthermore, (27) implies that the 
decay rate, which is equal to 7_, is independent of the initial amplitude a. This is 
shown in flgure [g] (c) where the time evolution of (Ag^(r)) is plotted for a = 1, 3/2 and 
3(72- = ^, 7-/72- = 1/10). 



4.3. Finite temperature^ NLD 



While at T = two-vibron losses facilitate the creation of the non-classical state (15), 



thermal two-vibron excitations will contribute to its decoherence. From (10) one sees 



that a flnite rate 72+ couples the states (n = 0, 1) to other excited states. These 
thermal excitations will eventually destroy the coherence of the non-classical states and 
the off-diagonal matrix elements will vanish. To quantify this decay we consider low 
temperatures, such that 72+ <C 72-. To obtain the relaxation rate to lowest order in 72+ 
it is suflicient to study excitations from n = to n = 2. Consequently, one gets from 



(10) the coupled equations 

?/^o(l, t) = (i/x - 472+)'0o + 272-V^2, 

7/^2(1, t) = 672+V^o + (5i// - I672+ - 472_)'02 



(28) 



For small 72+ we can 



supplemented by initial conditions for iPq{1,0) and ip2i^,0)- 
assume that these initial values correspond to the values of the state in (15), which is 
the stationary solution for 72+ = 0. This amounts to setting 

^o(l,t = 0) =^even(0), ^2(1,^ = 0) =0. (29) 



The system (28) can be solved for ipo{l,t), 

'372+ -i/U + 72- 



^o(l,t) 



cosh(2xt) + 
X exp[-(272- 



X 



sinh(2xt) 



(30) 



IO72+ - 3i/i)t]?/;even(0) , 



where X = [72- + (372+— 1/^)^+72- (972+— 2i/x)]-^/^. Equation (30) describes the relaxation 
of a non-classical state due to weak thermal excitations. The relaxation rate, which is 
determined by depends on the strength of the Kerr constant /x and temperature via 
72±. It is interesting to consider the limiting cases of strong and weak NLD. For strong 
NLD, 72-/m ~^ ^5 one flnds x ~^ 72- + (ll72+)/2 and ^00(1,^) becomes 

V;o(l,t)^V^o(l,0)e-^^+^ (31) 

Similarly, for weak NLD, 72-//^ ^ 0, x ^ — + 472+ + 72- and the time evolution is 

^o(l,t)^^o(l,0)e('^-^^^+)*. (32) 

Obviously, the decoherence rate increases with increasing strength of the Kerr constant. 
This can be understood from the fact that the accumulated phase shift in the excited 
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Figure 7. (a) Time dependence of the position variance of the nonhnearly damped 
mechanical mode at finite temperature for the same parameters as in figure |4] (b) 
Time dependence of |(Ag^(r)) — {Aq'^{oo))\T>o for parameters a = 1, /i = 72-, 7_ = 
and k^T/Vt G [1/4,1/2,3/4,1]. Symbols denote numerical results averaged over one 



period of and dashed lines show the behavior according to (35). (c) Time dependence 
of |/)o,i(r)P at kBT/Vt = 1/2 for a = 1 (o), a = 3/2 (□) and a = 3 (A). The long-time 
limit relaxation rate is according to (|34|. 



state is proportional to /x. Therefore, for finite /x there is an additional dephasing, which 
contributes to the decay of ipo{l^t). Altogether, the position variance becomes 

(Aq^) /ql ^]^+ 2n^{2Q) + Podd + |^o(l,t)r • (33) 

Figure [T] (a) shows the time dependence of the numerically calculated variance for 
the same parameters as in figure |3] (d) with the initial displacement amplitude a = 1. 
The T = variance is included as a reference (blue) and is almost identical to the 
curve for k^T /Vt = 1/4 (green). When comparing the graphs in figure [7] (a) with figure 
[5] (a) and |6] (a), one finds a qualitatively similar behavior. It is seen that also for 
finite temperatures the relaxation proceeds sequentially. The short-time dynamics show 
an increase in variance, corresponding to the formation of a Schrodinger cat state. 
For longer times the variance is increasing with time, and the rate of approaching 
the stationary value is faster for larger temperatures. In order to quantify the finite 
temperature relaxation rate, the time evolution of |(Ag^(r)) — (Ag^(oo))|T>o is shown 
in Figure [t] (b), where the stationary state variance (Ag^(oc))T>o is subtracted from 
(Ag^(r))T>o in order to remove the finite asymptotic value. The numerically obtained 
variance (symbols) is averaged over one period 27r///, which removes the oscillations due 
to the finite Kerr constant /x. The temperatures are the same as in figure [T] (a). The 
figure confirms the idea of different relaxation stages, with the decay in the second stage 
being exponential and the rate being temperature dependent. 

For comparison, figure [7|(c) shows the time evolution of |po, 1(^^)1^ for k^T/Vt = 1/2 
and a G [1,3/2,3]. It can be seen that in the long-time limit the decay rates do not 
depend on the initial oscillator displacement a. As discussed above, the off-diagonal 
matrix elements po,i decay due to thermal excitations. For low temperatures, their 
evolution is given by (30). When the Kerr constant and the NLD are equally strong, 
(/X = 72-), one finds x ~^ 72- + (1972+)/4 — i(/x — 372+/4) and the time evolution is 
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given by 

V;o(l, t) ^ 7/;o(l, o)ei(^+3^2+/2)^e-(^^^+/')^ (34) 



The resulting behavior is shown in figure [7| (c) and it can be seen that (34) correctly 
describes the decay of po,i- From figures [T] (a) and (b) it can be concluded that the 
variance increase is mainly due to the decay of |po,iP- Consequently, we make the 
Ansatz 



|(Ac?^(oo))t>o - (Ag^(r))|T>o = %'l^o(l, 0)| (35) 



which will be valid in the second stage of the relaxation. From (34) we find T{T) 



572+(2f]). This Ansatz is displayed in figure[7|(b) (dashed lines) and a good quantitative 



agreement with the numerical results is found for all temperatures. Using (35) one can 
therefore infer the value of |po,iP at the end of the first relaxation stage, which is 
dominated by the NLD. This allows for a reconstruction of the properties of the non- 



classical state (15), which is created at the end of the first relaxation stage. 



Finite temperature, LD and NLD 

Finally, in this section we consider the time evolution of the position variance at 
finite temperature under infiuence of both LD and NLD. Figure [s] (a) shows the time- 
dependence of the numerically obtained (Ag^) for 72- = /x, 7-/72- = 1/10, a = 1 and 
different temperatures. The behavior for T = 0, which corresponds to setting 7_ = 0, 
is included for reference (blue). As in the previous section we subtract the stationary 
state variance from (Ag^) to study the relaxation towards the stationary state. This 
is shown in figure M (b) . For long times an exponential decay is seen for temperatures 



k^T/Q < 1. To quantify the decay we use the Ansatz in (35) and consider 



\{Aq') - {Aq\oc))\T>o = qU^t - \)e-^\ (36) 

where a and F are estimated by a least-square fit to the numerical data. To this 
end the variance is averaged over one period 27r//i which removes the oscillations due 
to the finite Kerr constant /i. The fitting for T < 3/4 is done for the time interval 
1.5 • 10^ < r < 5.5 • 10^ and for T = 1 in the time interval of 10^ < r < 2 • 10^ The 
resulting behavior of ([36]) is displayed as dashed lines in figure [s] (b) . 

In figures [s] (c) and (d) the extracted parameters a and F are shown as functions of 
initial displacement amplitude a and temperature, respectively. Due to the sequential 
relaxation, the parameter a approximately corresponds to the variance of the non- 



classical state (15), which is created in the initial stage. The analytic expression 



for T = 0, given by equation (21), is shown in figure [8| (c) for comparison. The 
extracted low temperature results are in good agreement with the analytic curve. For 
higher temperatures deviations from the analytic behavior can be seen. However, the 



qualitative dependence on the initial amplitude is still similar to the one given by (21). 

In figure [s] (d) the temperature dependence of the extracted decay rates F is 
compared to the rates of one-vibron loss 7_ and two-vibron excitations 72+, which 
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Figure 8. (a) Time-dependence of the position variance of an initial coherent state 
with a = 1, at different temperatures with NLD (72- = /i) and LD (7_/72_ = 1/10). 
(b) Time dependence of |(Ag^(r)) — (Ag'^(oo))|T>o- The dashed Unes indicate a fit of 



(36) to the numerical data, (c) and (d) Extracted fit parameters a and F as functions 
of initial displacement amplitude a and temperature, respectively. The dashed line in 
(c) shows the variance according to (21) for T = and 7_ = 0. 



determine the decay rate in the situations investigated in sections |3.2| and |4.2[ It can be 
seen, that also in the present case the decay rate is independent of the initial coherent 
state displacement amplitude. For low temperatures, k^T/Vt < 0.3, one-vibron loss 
is the dominating damping mechanism. For higher temperatures the interplay of one- 
vibron and two-vibron processes leads to a more complicated behavior. At first the rate 
follows the behavior given by the combined rates 7_ + 672+ /2 but starts to deviate at 
k^T /Vt ^ 0.3. With increasing temperature it stays between 7_ +72+ and 7_ + 572+/2. 

The results presented in this section show that a ring-down measurement of the 
variance can be used to identify the quantum features of an initially created non-classical 
state, provided the sequential relaxation regime can be accessed. Moreover, the analysis 



Multi-phonon relaxation and generation of quantum states in a nonlinear mechanical oscillatorld 



of the decay rates can provide information about the dominating dissipation mechanisms 
in the system under investigation. 

5. Summary and conclusions 

In this work we have investigated a nonhnear osciUator in the quantum regime. The 
system is subject to nonhnear damping (two-vibron exchange), which resembles the 
two-photon decay known from quantum optics. The initiaUy displaced oscillator relaxes 
to a stationary state, which is determined by the relative strength of the linear and 
the nonlinear damping mechanisms and the temperature of the thermal bath. For 
sufficiently strong NLD we find that the relaxation is sequential: In the first stage. 



the initial density matrix decays to a non-classical state (15), followed by the second 
stage, where this state further decays to the final stationary state. The non-classical 
state depends on the initial state, while the decay rates depend on temperature and the 
strength of LD and NLD. 

At zero temperature, the sequential relaxation allows for a reconstruction of the 
properties of the non-classical state even in the presence of LD, i.e., when the stationary 
state is always the ground state. This can be achieved by performing a ring-down 
measurement of the position variance. For finite temperatures, NLD and thermal 
excitations lead to a different stationary state - a thermal state with a parity dependent 



Boltzmann-like distribution (23). The parity dependence is a signature of the NLD 
at finite temperatures. In the sequential relaxation regime the properties of the non- 
classical state can again be extracted from variance measurements. 

In the following we briefiy discuss the possibility of observing the relaxation of the 



non-classical state (15) in experiments with carbon-based nanomechanical resonators. 
A recent experiment |11| indicates significant NLD in micrometer sized, carbon-based 
resonators (nanotubes and graphene) with resonance frequencies of several hundreds 
of MHz. These systems are nonlinear and can be described by the classical Duffing 
equation of motion [36j 

mq —mVt^q — a^q^ — iq — Vq'^Q^ (37) 

where q is the oscillator position, m and Q are the oscillator mass and frequency, ao 
is the nonlinearity strength, 7 is the LD rate and rj is the NLD rate. By comparing 



(37) with the corresponding quantum mechanical equation, the relation between the 



parameters (/x,7_,72_) and the parameters of the Duffing equation can be established, 
^ 3hao ^ hr] ^ ^ 

The temperature for reaching the quantum regime without cooling is of the order of 
T < 10 mK for resonators with frequencies in the 100 MHz range. The values obtained 
correspond to /x/f] <C 10~^, meaning that the RWA-analysis works well. The 



m 



11 



magnitude of the nonlinear damping rate is found to yield 72- ^ /x, which is promising 
for the formation of a non-classical state. The condition to find the sequential relaxation 
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regime amounts to 

7- < 72- < — > 7 < 



hr] 



1 



(39) 



where qo is the zero point amphtude. Consequently, for systems with a large zero 
point amplitude it is easier to reach the sequential relaxation regime. In this regard, 
carbon-based resonators are very promising due to their low mass and tunable frequency. 
In present experiments the condition 7_ < 72- is not yet fulfilled, but the ongoing 
experimental and technological progress in the field will help to improve the parameters 



in order to meet the condition (39). 



In conclusion, our study shows that exploiting NLD present in (carbon-based) 
NEMS resonators is a promising route to generation and investigation of non-classical 
effects in NEMS. Studies of the sequential relaxation in multi-partite systems, such 
as coupled oscillators, would be of interest with respect to the generation of quantum 
entanglement. 
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Appendix A. Derivation of equations (17) and (19) 



For zero temperature we set 72+ = and get 
d 

Q^PnAt) = l2-{n + 2)(n + l)Pn+2,n+2 " 72-^(^ " ^) Pn 

d 



(A.l) 



dt 



t) = 72- (n + 2)(n + l)V^,+2 - + 2n) + 72_n'] • (A.2) 



One sees that for t ^ oc only p^o, and ^00(1) remain non-zero (depending on initial 
conditions). To find the respective values in the stationary limit, one first verifies that 



d ^ 



2n 







d 



n+l,2n+l 







(A.3) 



Therefore, Po,o(oo) = Z;„p2n,2n(0) and pi,i(oo) = P2n+i,2n+i(0) are constants of 
motion. There exists another constant of the motion, which can be found as follows: 
We take tpn = e^'^^tpn and consider 

d ( 

J] C2n^2n = 72- + 2)(2n + l)i'2n+2 " 



dt 



n=0 



-i— 4n + 4n^ 

72- 



i'2n (A. 



with so far unknown coefficients C2n- Those shall be determined such that the right 
hand side vanishes. This is the case if 



C2n-2(2n-l) 



u 

-12— + 2n 

72- 



C, 



2n 



(A.5) 
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is fulfilled. This recursion can be solved in terms of Gamma functions 

r(| + n)r(i-i^) 

0rr(l + n-i:^) 

which can be verified using the identity T{z + 1) = zT{z). Hence, ■(/'o(l,oo) = 

EnC2n^2n(0). 

If the initial state is a coherent state, 



pn,m{0) = exp(-|a|2)a"a*"^/Vn!m! , (A.7) 

one gets 



J]C2n^2n(0) = exp(-|a|2)a* J]C2n7^ =exp(-|a|2)a*F(;^) (A.8) 

with z = It can be shown, that F{z) corresponds to the definition of the confiuent 
hypergeometric hmit function qFi [37]. Consequently, 

Ml, oo) = V C2nM{0) = exp{-\a\')a*oF,{l - i-^, |a|V4) . (A.9) 

^ 72- 



Appendix B. Derivation of equations (23) and (24) 



The equations for finite temperature can be directly obtained from (jTol): 



§-tPn,n{t) = 'y2+n{n - l)pn-2,n-2+72-{n + 2){n+l)pn+2,n+2 

- [72-n(n - 1) + 72+{n + l)(n + 2)] pn,n , 
f^ipn{l,t) = 72+n(n + l)t/;„_2 + 72-(n + 2)(n + l)t/;„+2 

- [-i^(l + 2n) + 72-n2 ^ ^^^(^ ^ 2)2] . 



(B.l) 



The stationary probabihty distribution Pn,n(co) obeys the detailed balance condition 

72+Pn,n(oc) = 72-Pn+2,n+2 (oc) . (B.2) 

Therefore, 

P2n,2n(00) = (^^^ po,o(oo) = e'^^^'^Vo.O (oo) , (B.3) 

\ n 

72+ \ / X ^-2/if7/3n^ 



P2n+l,2n+l(00) = j pi,i(oo) = e'^'^'^'^Vi,! (oo) . 

Since parity is conserved one can use 

P2n,2n(00) = Po,o(^) J] e-^^*^^" = Peven , (B.4) 

n n 

Y,P2n+l,2n+l{oc) = (^) ^ e'^'^^^" = Podd 

n n 

to fix the values po,o and pi^i. The sum evaluates to ^^e~^^^^^ = 1 + n5(2f]). With 
these results it is straight-forward to calculate the average number of excitations 
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^^''"^^ -nB{2flf ( + 1 ) = 2PevennB(21^) 



l + nB{2n) ' \nB{2n) 

2 P 

V(2n + l)p2n+i,2n+i(oc) = ^-^^ V ne"^^^^- + Podd = PoM{2nB{2n) + 1). 

And therefore 

^a'^a) = y]2np2n,2n(oQ)+y](2n+l)p2n+i,2n+i(oQ) = 2nB(2f^)+Podd .(B.5) 

2n 2n 
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